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I.   ;ntroduction  and  Description  of  Tests  Performed 

The  purpose  of  this  report  is  to  present  the  results  of  various 
tests  made  concerning  the  numerical  integration  of  differential  equations. 
The  routine  used  in  the  tests  is  a  Fortran  program  which  integrates  a  set 
of  N  first  order  ordinary  differential  equations  using  predictor-corrector 
multistep  methods.   The  theoretical  basis  for  the  methods  used  is  given  in 
Report  No.  221,  Numerical  Integration  of  Stiff  Ordinary  Differential 
Equations ,  by  C.  W.  Gear.   The  coefficients  for  these  methods  are  given  in 
Appendix  I. 

The  main  program  uses  two  subroutines  which  are  related  specifically 
to  the  equations  being  integrated.   The  first  evaluates  the  function  which 
is  given  in  the  following  form: 
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Where  h  is  the  current  step  size.   The  second  of  these  subroutines 
calculates  the  following  matrix: 
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Where  a  is  the  leading  coefficient  for  the  method  being  used.   This 
matrix  is  used  in  obtaining  the  corrections  to  the  predicted  values  of 
the  function. 


The  initial  values  of  the  function  must  be  given,  along  with  the 
upper  bound  on  the  error  (absolute  error  is  used),  the  number  of  corrector 
iterations  to  be  allowed  in  one  attempt  to  integrate,  and  the  original  step 
size  to  be  used.   The  step  size  is  thereafter  automatically  controlled. 
The  criteria  for  convergence  and  doubling  or  halving  of  the  step  size  are 
derived  from  the  error  bound  provided  initially. 

The  tests  involved  the  use  of  two  versions  of  the  integration 
program,  one  in  which  the  order  of  the  method  to  be  used  is  specified,  and 
another  in  which  the  order  is  automatically  chosen  and  changed  depending  on 
the  course  of  the  integration.   The  order  of  the  method  is  subject  to 
change  if  the  step  size  is  halved  or  doubled.   The  decision  to  change  order 
or  not  is  based  on  the  comparison  of  estimates  of  the  error  in  various 
order  methods . 

The  purpose  of  the  tests  which  were  made  was  as  follows: 

1.  The  investigation  of  the  efficiency  of  various 
numbers  of  corrector  iterations  allowed  with 
various  order  methods. 

2.  The  comparison  of  stiff  and  non-stiff  methods 
on  both  (stiff  and  non-stiff)  types  of 
equations  in  terms  of  accuracy  and  time  cost. 

3.  The  comparison  of  various  fixed  order  methods 
and  automatic  order  selection  with  respect  to 
cost  and  accuracy. 

The  description  of  the  four  systems  of  equations  used  for  the 
tests  may  be  found  in  Appendix  II. 


II.   The  Question  of  the  Number  of  Correction  Iterations 

It  can  be  seen  from  the  graph  in  Appendix  III  that  it  is  in 
general  not  advisable  to  use  four  corrector  iterations,  as  this  is  as 
costly  or  more  costly  than  using  two  or  three  iterations.   The  decision 


whether  to  use  two  or  three  iterations  is  a  rather  difficult  one  to  make. 
The  results  of  this  comparison  made  on  the  other  equation  sets  are 
essentially  the  same  as  equation  set  I  (shown  in  the  graph),  though  the 
use  of  three  iterations  with  the  other  equations  seems  to  be  more 
advantageous  in  some  cases.   The  only  marked  differences  between  two  and 
three  iterations  occurs  in  the  high  order  method  where  three  iterations 
i  s  preferable  for  large  and  medium  errors  and  two  iterations  for  small 
errors.   The  advantage  in  using  two  iterations  in  the  high  order-small 
error  case  however  was  not  indicated  in  tests  on  the  other  equation  sets. 
In  general  two  iterations  is  more  efficient  for  low  order-small  error 
and  three  iterations  is  better  for  high  order-large  error  situations. 

The  value  of  using  two  or  three  iterations  seems  to  depend  on 
the  equations  being  integrated  to  some  extent,  but  on  the  basis  of  these 
results,  the  remaining  tests  were  performed  using  three  corrector 
iterations  exclusively. 

A  remark  is  in  order  concerning  the  errors  indicated  on  this 
graph.   These  are  not  the  actual  errors,  but  the  error  bound  provided 
for  the  integration  program.   The  actual  error  and  desired  error  correspond 
closely  for  the  cases  where;  l)  low  order  is  used  with  a  large  error  bound, 
2)  medium  order  is  used  with  a  medium  error  bound,  3)  and  high  order  is 
used  with  a  small  error  bound.   When  low  order  methods  were  used  with  a 
small  error  bound,  the  actual  error  exceeded  the  bound  given.   However 
when  high  order  methods  were  used  with  a  large  error  bound,  more  accuracy 
was  obtained  than  was  requested. 

It  can  be  seen  that  the  two  iteration  process  uses  generally  the 
fewest  number  of  function  evaluations  and  the  greatest  number  of  matrix 
evaluations,  while  four  iterations  uses  the  greatest  number  of  function 
evaluations  and  the  fewest  matrix  evaluations.   The  use  of  three  iterations 
usually  falls  between  the  other  two  but  sometimes  requires  the  fewest 
function  evaluations  or  matrix  evaluations . 

It  was  also  seen  that  the  number  of  corrector  iterations  used  is 
not  consistently  related  to  the  actual  error  produced.   That  is,  three 


or  four  corrector  iterations  would  not  in  general  yield  more  accurate 
results  than  two  iterations  for  example. 

III.   Stiff  Versus  Non-Stiff  Integration  Methods 

The  four  graphs  in  Appendix  IV  show,  for  each  equation  set,  the 
relationship  between  the  actual  error  and  the  cost  of  the  integration  in 
terms  of  function  evaluations  required.   The  results  are  shown  for  stiff 
and  non-stiff  methods  of  order  two,  four,  and  six,  and  automatic  order 
selection.   Three  corrector  iterations  were  allowed  in  every  case. 

In  general  it  can  be  remarked  that  the  graph  of  the  method  of 
order  two  has  the  greatest  slope,  while  the  graph  of  the  sixth  order 
method  has  the  smallest.   That  is,  it  is  preferable  to  use  low  order 
methods  when  a  large  error  can  be  tolerated,  and  high  order  methods  when 
more  accuracy  is  desired. 

The  graphs  of  the  stiffly  stable  equations  show  that  non-stiff 
methods  used  on  stiffly  stable  equations  are  in  all  cases  inferior  to  the 
stiff  methods.   The  integration  was  not  actually  performed  over  the  entire 
intervals  of  integration  indicated  in  the  description  of  the  equation  sets, 
since  it  was  found  to  be  very  impractical  in  terms  of  computer  time 
expended.   The  results  shown  are  an  extrapolation  based  on  the  assumption 
that  the  actual  error  would  not  have  been  substantially  larger  if  the 
integration  interval  had  been  completely  traversed. 

The  graphs  of  the  non-stiff  equations  however,  show  that  the 
non-stiff  methods  are  superior  to  the  stiff  methods  of  the  same  order, 
but  that  some  stiff  methods  are  better  than  the  non-stiff  methods  of  a 
different  order.   It  is  also  evident  that  in  most  cases  the  differences 
in  the  two  types  of  methods  for  the  same  order  are  not  appreciable.   It 
is  therefore  reasonable  to  choose  stiffly  stable  methods  over  non-stiff 
methods  when  integrating  an  arbitarary  system  of  equations  which  may  be 
stiff  or  non-stiff. 


IV.   Fixed  and  Varying  Order  Integration  Techniques 

In  making  the  comparison  between  fixed  order  and  automatic  order 
selection,  it  will  be  helpful  to  discuss  the  stiff  and  non-stiff  equation 
sets  separately,  since  the  behavior  of  each  group  is  dependent  on  the 
coefficients  (stiff  or  non-stiff),  used  with  the  automatic  order  selection 
program. 

Equation  sets  I  and  II  (stiffly  stable)  were  integrated  more 
successfully  in  general  using  the  stiff  automatic  order  selection  version 
than  when  the  fixed  order  methods  were  used,  considering  the  range  of 
error  as  a  whole.   The  exception  to  this  is  that  stiff  fixed  order  methods 
were  sometimes  more  effective  when  a  relatively  large  error  was  allowed. 

With  equation  set  I,  the  non-stiff  automatic  order  selection 
version  was  only  slightly  better  than  the  fixed  order  six  method,  and 
was  more  costly  than  all  other  methods. 

The  behavior  of  equation  set  II  is  more  complicated.   The  last 
remarks  concerning  equation  set  I  apply  to  equation  set  II  only  in  the 
small  error  range.   For  medium-sized  errors,  the  non-stiff  automatic 
order  selection  method  is  more  effective  than  non-stiff  orders  four  and 
six,  and  slightly  less  effective  than  order  two.   The  result  is  a  plateau- 
like graph.   There  is  another  sharp  drop  (not  shown  on  the  graph), 
corresponding  to  the  transition  from  medium  to  large  error.   The  explanation 
for  this  phenomenon  seems  to  be  that  the  high  orders  (five  and  six)  are 
chosen  by  the  order-selecting  algorithm  when  small  error  is  requested, 
while  lower  orders  predominate  when  medium  or  large  error  is  permitted. 

The  integration  of  the  non-stiff  equations  indicates  that  automatic 
order  selection  (using  stiff  methods)  is  more  effective  than  the  fixed 
order  process  using  stiff  methods,  and  approximately  equivalent  to  or 
better  than  the  non-stiff  methods.   The  non-stiff  automatic  order  selection 
program  is  superior  to  the  stiff  automatic  version  for  both  equation  sets 
III  and  IV,  especially  in  the  large  error  range.   However  the  non-stiff 
automatic  method  is  in  general  not  quite  as  good  as  the  fixed  order  six 
method. 
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These  results  would  suggest  that  for  the  integration  of  a  general 
system  of  ordinary  differential  equations  where  it  is  not  known  whether 
the  system  is  stiffly  stable  or  not,  it  is  more  efficient  to  use  stiffly 
stable  integration  methods.   Also  it  is  indicated  that  it  is  preferable 
to  vary  the  order  of  the  method  being  employed  during  the  course  of  the 
integration,  especially  when  using  stiff  methods. 

V.   Some  Additional  Results 


Since  the  tests  described  thus  far  have  been  completed,  the 
automatic  order  selection  integration  program  has  been  revised.   The 
significant  changes  and  a  comparison  between  the  results  using  the  old 
and  new  version  are  found  in  Appendix  V. 

Results  are  shown  for  the  new  version  for  two  different  cases. 
The  first  case  uses  a  constant ,10  ,  as  the  criterion  for  convergence 
•of  the  corrector.   The  second  case  uses  the  variable 


(-P-3) 
E  =  2       *  requested  error 


where  P  is  the  order  of  the  method  being  used.   In  the  first  case  it  can 
be  seen  that  there  is  little  correspondence  between  the  requested  and 
actual  errors,  and  much  more  accuracy  is  achieved  than  is  requested  in 
many  instances.   This  is  true  of  course  because  the  bound  does  not  depend 
on  the  requested  error. 

When  the  new  version  was  used  with  the  variable  criterion  for 
convergence  (which  is  the  same  as  that  used  with  the  old  version)  the 
correspondence  between  requested  and  actual  error  is  essentially  restored. 
This  would  seem  to  be  a  desirable  feature,  even  though  in  a  few  cases  the 
constant  convergence  criterion  seems  to  be  more  economical. 

Overall,  it  can  be  seen  that  the  present  program  is  less  costly 
then  the  previous  version  for  a  comparable  amount  of  accuracy. 


The  partial  derivative  matrix  has,  in  all  of  these  tests, 
been  evaluated  using  the  values  of  the  solution  at  the  last  successful 
step  in  the  integration.   It  has  been  found  that  a  decrease  in  the  number 
of  function  evaluations  necessary  can  be  achieved  by  evaluating  the  matrix 
using  the  predicted  values  of  the  solution.   Appendix  VI  contains  a 
comparison  of  these  two  alternatives. 


Stiff  Methods 


APPENDIX  I 
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APPENDIX  II 

DESCRIPTION  OF  SYSTEMS  OF  EQUATIONS  TESTED 

y^  =  ((R  -  Ay3  -E^  +  Ey2)  /  V  y^Q)  =  1 

J2   =  W(yi  -  y2)  y2(0)  =  1 


y^  =  Pyx  /  C  y3(0)  =  0 


A  =  .0001  E  =  .0065  V  =  .0001 

C  =  200.  P  =  20.  W  =  .0785 


This  system  is  stiffly  stable  and  was  integrated  from 
x  =  0  to  x  =  500. 
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II-      y-L  =  -  Ay-L  -  Byxy3  y^O)  =  1 


y2  =  -  Cygy  y2(0)  =  1 


y3  =  Ayx  -  Byxy3  -  Cy2y3  y^O)  =  0 


A  =  .013 
B  =  1000. 

C  =  2500. 

This  set  is  also  stiffly  stable  and  integrated  from 
x  =  0  to  x  =  50. 


III.     y  =  y2y  y1(0)  =  0 


y2  -  -  y-^  y2(°)  =  i 


y3  =  -  .517^2  y3(0)  =  1 


This  system  is  not  stiffly  stable  and  was  integrated 
from  x  =  0  to  x  =  10. 


IV. 


y-L  =  v1  yx(o)  =  i 


y2  =  -  y2  y2(o)  =  i 


y3  =  yk  y3(o)  =  o 


y^  =  -  y3  yu(o)  =  l 

This  system  is  not  stiffly  stable  and  was  also  integrated 
from  x  =  0  to  x  =  10. 
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APPENDIX  V 

LIST  OF  SIGNIFICANT  CHANGES  IN  THE  AUTOMATIC  ORDER 
SELECTION  INTEGRATION  PROGRAM 


Old  Version 


New  Version 


The  step  size  can  only  be  changed 
by  a  factor  of  2  or  .5  at  any  one 
time. 


The  step  size  can  be  changed  by  a 
factor  of  2n  (n  is  any  positive  or 
negative  integer).   This  is  limited 
only  by  a  control  involving  the 
printing  points  for  results. 


Old  error  coefficients  are  used 
in  order-changing  criteria  only. 


New  error  coefficients  are  used  for 
order-changing  criteria  and  for  the 
test  for  successful  integration. 


The  test  for  successful  integration 
(error  bound)  is  simply  the 
requested  error. 


The  error  bound  on  the  integration 
depends  on  requested  error,  error 
coefficients,  the  order  of  method 
being  used,  and  step  size. 


Error  bound  for  integration  is  used 
as  follows : 


Error  bound  is  used  as  follows 


Error  Bound  > 


predicted   corrected 
yi        -yi 


n( Error  Bound)  >L 
i=l 


P    c 
y.  -  y. 


for  all  i  in  order  for  step  to  be 
taken. 


where  n  is  the  number  of  equations 


Order  can  be  changed  after  step 
doubling,  halving  (non-convergence), 
partial  derivative  matrix  evaluation, 
and  it  may  be  changed  upward  or 
downward. 


Order  can  be  changed  only  after  a 
certain  number  of  steps  (depending 
on  the  order),  and  if  non-convergent, 
the  order  can  only  be  decreased  if 
changed. 


If  the  corrector  has  not  converged, 
step  halving  or  matrix  re-evaluation 
is  performed  in  the  main  program  and 
order  changing  is  allowed  immediately. 


If  the  corrector  does  not  converge, 
alternate  matrix  re-evaluation  and 
step  size  changing  (by  a  factor  of 
2~j)  take  place  in  the  integration 
subroutine',  until  convergence  is 
accomplished  or  the  step  size  becomes 
too  small. 


APPENDIX  VI 

COMPARISON  OF  MATRIX  EVALUATION  USING  LAST  ACCEPTED  VALUES 
OF  THE  SOLUTION  AND, USING  PREDICTED  VALUES  OF  THE  SOLUTION 
(ACTUAL  ERROR/0  OF  FN.  EVALUATIONS) 


REQ. 
ERROR 

EQUATION 
SET   I 

EQUATION 
SET  II 

1 

EQUATION 

SETni 

EQUATION 
SET  IV 

A 
ID"3 

p 

i 
7.5xlO-2 
315 

2.1+xlO-3 
Ikk 

UxlO~3 
260 

l.UxlO° 

27h 

7.1xl0"2 
308 

2.7xlO_3 
1*5 

2.6xl0~3 
219 

l.UxlO0 
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A 
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P 
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620 

3.8x10  p 
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551   ■ 
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309 

Q 
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986 
l 

,   3-h-xIO"5 
1191 

=  Accepted 
=  Predicted 


COMPARISON  OF  OLD  AND  NEW  VERSIONS  OF  THE  AUTOMATIC  ORDER 
SELECTION  PROGRAM  IN  TERMS  OF: 
(ACTUAL  ERROR/ #  OF  FN.  EVALUATIONS) 
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